#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define UNROLL_FACTOR 4  // 循环展开因子

void matrix_multiply_unrolled(double *A, double *B, double *C, int m, int n, int k) {
    for (int i = 0; i < m; i++) {
        for (int l = 0; l < n; l++) {
            double a = A[i * n + l];  // 将A的元素存入寄存器减少访问次数
            
            // 主循环（展开部分）
            int j = 0;
            for (; j <= k - UNROLL_FACTOR; j += UNROLL_FACTOR) {
                C[i * k + j]     += a * B[l * k + j];
                C[i * k + j + 1] += a * B[l * k + j + 1];
                C[i * k + j + 2] += a * B[l * k + j + 2];
                C[i * k + j + 3] += a * B[l * k + j + 3];
            }
            
            // 处理剩余部分（尾循环）
            for (; j < k; j++) {
                C[i * k + j] += a * B[l * k + j];
            }
        }
    }
}

int main() {
    const int m = 512, n = 1024, k = 2048;
    
    // 使用对齐内存分配（可选）
    double *A = (double *)aligned_alloc(64, m * n * sizeof(double));
    double *B = (double *)aligned_alloc(64, n * k * sizeof(double));
    double *C = (double *)aligned_alloc(64, m * k * sizeof(double));

    // 初始化随机数生成器
    srand(time(NULL));
    
    // 初始化矩阵（使用循环展开优化初始化）
    for (int i = 0; i < m * n; i++) A[i] = (double)rand() / RAND_MAX;
    for (int i = 0; i < n * k; i++) B[i] = (double)rand() / RAND_MAX;
    for (int i = 0; i < m * k; i++) C[i] = 0.0;

    // 计时
    clock_t start = clock();
    matrix_multiply_unrolled(A, B, C, m, n, k);
    clock_t end = clock();

    double time_taken = ((double)(end - start)) / CLOCKS_PER_SEC;
    printf("Optimized C Time: %.6f seconds\n", time_taken);

    // 计算FLOPS
    double flops = 2.0 * m * n * k / time_taken / 1e9;
    printf("Performance: %.2f GFLOPS\n", flops);

    free(A);
    free(B);
    free(C);
    return 0;
}